#lang racket

(require math/bigfloat)

(define epsilon (bf 0.0000000001))
(define max-denom (bf/ 1.bf epsilon))
(define cf-bound 100)
(define n-bound 100)

(define (bflogxy x y) (bf/ (bflog y) (bflog x)))
(define (bfpowxy x y) (bfexp (bf* y (bflog x))))
(define (bfave x y) (bf/ (bf+ x y) 2.bf))

(define (continued-fraction bfinput n)
  (cf-iter bfinput n n))

(define (cf-iter c i n)
  (if (bf> (bfabs c) epsilon)
      (let* ((a (bffloor c))
             (nextc (bf/ 1.bf (bf- c a))))
        (if (or (bf> nextc max-denom) (= i 0)) (cons a '())
            (cons a (cf-iter nextc (- i 1) n))))
      '()))

(define (cf-calc ak cfk-1 cfk-2)
  (define pk-1 (car cfk-1))
  (define qk-1 (cdr cfk-1))
  (define pk-2 (car cfk-2))
  (define qk-2 (cdr cfk-2))
  (define pk (bf+ (bf* ak pk-1) pk-2))
  (define qk (bf+ (bf* ak qk-1) qk-2))
  (cons pk qk))

(define (cf-approxy cf k)
  (define n (min k (length cf)))
  (cond
    ((= n 0) (cons (car cf) 1.bf))
    ((= n 1) (cons (bf+ (bf* (car cf) (cadr cf)) 1.bf) (cadr cf)))
    (else (cf-calc (list-ref cf n)
                   (cf-approxy cf (- n 1))
                   (cf-approxy cf (- n 2))))))

;(continued-fraction pi.bf 10)
;(cf-approxy (continued-fraction pi.bf 10) 4)

(define (factoried s t N n)
  (define N4 (bf* s t N 4.bf))
  (define sqrtN4 (bfsqrt N4))
  (define M (bfceiling sqrtN4))
  (define M+n (bf+ M n))
  (define M2 (bf* M+n M+n))
  (define delta (bf- M2 N4))
  (define sqrtdelta (bfsqrt delta))
  (if (bf= n 0.bf) #f
      (if (bf= sqrtdelta (bffloor sqrtdelta))
          (cons (bf/ (bfave M+n sqrtdelta) t) (bf/ (bfave M+n (bf- sqrtdelta)) s))
          (factoried s t N (bf- n 1.bf)))))

(define (factoring-recur N n t)
  (define sqrtN (bfsqrt N))
  (displayln sqrtN) ;ok
  (define qpratio (bfpowxy N (bf/ t (bfpowxy 2.bf n))))
  (displayln qpratio) ;ok
  (define alpha (bf/ (bf+ qpratio 1.bf) (bf- qpratio 1.bf)))
  (displayln alpha) ;ok
  (define next-qpratio (bfsqrt qpratio))
  (displayln next-qpratio) ;ok
  (define next-alpha (bf/ (bf+ next-qpratio 1.bf) (bf- next-qpratio 1.bf)))
  (displayln next-alpha) ;ok
  (define beta (bfave alpha next-alpha))
  (displayln beta) ;ok
  (define cf-beta (continued-fraction beta cf-bound))
  ;(displayln cf-beta) ;ok
  (define cf-approxy0 (cons (car cf-beta) 1.bf))
  (define cf-approxy1 (cons (bf+ (bf* (car cf-beta) (cadr cf-beta)) 1.bf) (cadr cf-beta)))
  (define cf-tail (cddr cf-beta))
  ;(displayln cf-approxy0) ;ok
  ;(displayln cf-approxy1) ;ok
  ;(displayln cf-tail) ;ok
  ;(define cfk-2 cf-approxy0)
  ;(define cfk-1 cf-approxy1)
;  (define (cf-approxy-increase)
;    (if (null? cf-tail) '()
;        (let ((cf-next (cf-calc (car cf-tail) cfk-1 cfk-2)))
;          (set! cfk-2 cfk-1)
;          (set! cfk-1 cf-next)
;          (set! cf-tail (cdr cf-tail))))
;    cfk-1)

  (define cf-current (list cf-tail cf-approxy1 cf-approxy0))

  (define (cf-approxy-increase cf-cur)
    (if (null? cf-tail) '()
        (let* ((ak (first cf-cur))
               (cfk-1 (second cf-cur))
               (cfk-2 (third cf-cur))
               (cfk (cf-calc (car ak) cfk-1 cfk-2)))
          (list (cdr ak) cfk cfk-1))))

  

  ;(displayln (cf-approxy-increase cfcur)) ;ok
  ;(displayln (cf-approxy-increase (cf-approxy-increase cfcur))) ;ok
  ;(displayln (cf-approxy-increase (cf-approxy-increase (cf-approxy-increase cfcur)))) ;ok

  (define (cf-approxy-sqrtN cf-cur denom)
    (define cf-tail (first cf-cur))
    (define cfk (second cf-cur))
    (if (and (bf< (cdr cfk) denom) (not (null? cf-tail)))
        (let ((cf-next (cf-approxy-increase cf-cur)))
          (cf-approxy-sqrtN cf-next denom))
        cf-cur))

  (define cf-approxy (second (cf-approxy-increase (cf-approxy-sqrtN cf-current sqrtN))))

  (displayln cf-approxy) ;ok

  (define (f a c d n-try)
    (define b (bf/ c d))
    (define b-a (bf- b a))
    (define a2-1 (bf- (bf* a a) 1.bf))
    (define ab-1 (bf- (bf* a b) 1.bf))
    (define rho (bf/ (bf* n-try (bfsqrt a2-1)) (bf* d sqrtN)))
    ;(bf- (bf* b-a b-a) (bf* rho ab-1) (bf/ (bf* rho rho) 4.bf))
    (bf+ (bf* rho ab-1) (bf/ (bf* rho rho) 4.bf))

    )

  (displayln (f (bf 5.8) (bf 45718) (bf 10357) (bf 100)))
  
  #|
  ; (define beta-approxy (bf/ (car cf-approxy) (cdr cf-approxy)))
  (define (f a c d n-try)
    (define b (bf/ c d))
    (define b-a (bf- b a))
    (define a2-1 (bf- (bf* a a) 1.bf))
    (define ab-1 (bf- (bf* a b) 1.bf))
    (define rho (bf/ (bf* n-try (bfsqrt a2-1)) (bf* d sqrtN)))
    (bf- (bf* b-a b-a) (bf* rho ab-1) (bf/ (bf* rho rho) 4.bf)))
  (define (find-n i cf)
    (define c (car cf))
    (define d (cdr cf))

;    (displayln i)
;    (displayln cf)
;    (displayln (f alpha c d i))
;    (displayln (f next-alpha c d i))

    
    (if (and (bf< (f alpha c d i) 0.bf) (bf< (f next-alpha c d i) 0.bf)) (bf i)
        (if (bf< i (bf n-bound)) (find-n (bf+ i 1.bf) cf)
            0.bf)))
  (define (find-beta cfk)
    (define n-try (find-n 1.bf cfk))
    (if (bf= n-try 0.bf)
        (find-beta (cf-approxy-increase))
        n-try))
  ;(define n-try (find-beta cf-approxy))

  (define c (car cfk-1))
  (define d (cdr cfk-1))

;  (displayln c)
;  (displayln d)
;  (displayln n-try)
  ;(define r (factoried (bf+ c d) (bf- c d) N (bf n-try)))

;  (if r (displayln r)
;      (factoring-recur N n (bf/ t 2.bf)))

  |#

  '()
  )

(define (factoring1 N)
  (define n (bfceiling (bf- (bflog2 (bf- 1.bf (bflogxy N (bf- N 1.bf)))))))
  (displayln n)
  (define t (bfceiling (bf/ (bfpowxy 2.bf n) (bflog2 N))))
  (displayln t)
  (factoring-recur N n t))

;(factoring (bf* (bf 1087) (bf 1493)))


(define (factoring N n)
  (define sqrtN (bfsqrt N))
  (define beta-inc (bf/ pi.bf (bf 360)))
  (define beta (bf* (bf/ n (bflog n)) beta-inc))
  (define cos-beta (bfcos beta))
  (define cf-beta (continued-fraction cos-beta cf-bound))
  (define cf-approxy0 (cons (car cf-beta) 1.bf))
  (define cf-approxy1 (cons (bf+ (bf* (car cf-beta) (cadr cf-beta)) 1.bf) (cadr cf-beta)))
  (define cf-tail (cddr cf-beta))
  (define cf-current (list cf-tail cf-approxy1 cf-approxy0))

  (define (cf-approxy-increase cf-cur)
    (if (null? (first cf-cur)) cf-current
        (let* ((ak (first cf-cur))
               (cfk-1 (second cf-cur))
               (cfk-2 (third cf-cur))
               (cfk (cf-calc (car ak) cfk-1 cfk-2)))
          (list (cdr ak) cfk cfk-1))))
  (define (cf-approxy-sqrtN cf-cur denom)
    (define cf-tail (first cf-cur))
    (define cfk (second cf-cur))
    (if (and (bf< (cdr cfk) denom) (not (null? cf-tail)))
        (let ((cf-next (cf-approxy-increase cf-cur)))
          (cf-approxy-sqrtN cf-next denom))
        cf-cur))

  (define cf-approxy (second (cf-approxy-increase (cf-approxy-sqrtN cf-current sqrtN))))

  (define c (car cf-approxy))
  (define d (cdr cf-approxy))
;  (displayln (bf/ 1.bf cos-beta))
;  (displayln c)
;  (displayln d)

  (define r (factoried (bf+ c d) (bf- c d) N (bf 100)))
  (displayln r)
  
  (if (bf< n (bf* (bf 179) (bflog (bf 180)))) (factoring N (bf+ n 1.bf))
      '()))

;(factoring (bf* (bf 1087) (bf 1493)) (bf 2))

(define (bfiterate f start end)
  (if (bf<= start end)
      (begin
        (display start)
        (display #\space)
        (displayln (f start))
        (bfiterate f (bf+ start (bf 1.0)) end))
      '()))

(define (gougu n)
  (define 2n+1 (bf+ (bf* 2.bf n) 1.bf))
  (bf/ (bf+ (bf* 2n+1 2n+1) 1.bf) 2.bf))

;(bfiterate gougu 1.bf (bf 1000))

(define (fenjie n)
  (define sqrtn (bfsqrt n))
  (define a (bfceiling sqrtn))
  (define theta (bf- a sqrtn))
  (define delta (bf- (bf* a a) n))
  (define sqrt-delta (bffloor (bfsqrt delta)))
  (if (bf= delta (bf* sqrt-delta sqrt-delta))
      (list a sqrt-delta theta)
      #f))

;(cf-approxy (continued-fraction (bf 1.0143360285728458752914331068009) 10) 4)

;(bfiterate fenjie 1.bf (bf 10000))

(define (search N)
  (define theN (bf* N (bf 20159)))
  (define sqrtN (bfsqrt theN))
  (define M (bfceiling sqrtN))
  (define theta (bf- M sqrtN))
  (if (and (bf> theta (bf/ 1.bf (bfsqrt sqrtN)))
           (bf< (bf- (bf* M M) theN) (bf- (bf* M 2.bf) 1.bf)))
      N
      '()))


(bfiterate search (bf 1) (bf 1000))



